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ABSTRACT 


Numerical  solutions  of  the  steady-state  Navier-Stokes  equations  are 
obtained  which  describe  flow  past  a  square  protuberance  immersed  in  a 
plane  Couette  flow.  The  numerical  solutions  are  illustrated  with  velo¬ 
city  vector  plots,  and  plots  of  streamlines  and  vorticity  contours  in 
the  flow  field  for  Reynolds  numbers  between  1  and  200,  based  on  plate 
velocity  and  protuberance  height.  A  small  separation  bubble  is  oredicted 
upstream  of  the  protuberance  with  length  and  height  almost  independent 
of  Reynolds  number.  The  downstream  separation  bubble  has  a  height 
always  less  than  the  protuberance  height,  but  a  length  which  increases 
almost  linearly  with  Reynolds  number  reaching  more  than  twenty-five 
protuberance  heights  at  Re  =  200.  The  flow  perturbations  caused  by  the 
protuberance  are  shown  to  persist  far  downstream  of  the  protuberance  at 
Re  =  100  and  200. 
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I,  INTRODUCTION 


The  qualitative  features  of  the  flow  field  downstream  of  two- 
dimensional  roughness  elements  immersed  in  a  laminar  flat  nlate  boundary 
layer  are  well  known.  In  many  wind  tunnel  experiments,  such  protuberances 
are  employed  to  promote  transition  and  assure  turbulent  boundary  layer 
characteristics  on  models.  At  high  Reynolds  number,  the  flow  at  the 
base  of  the  protuberance  is  separated  with  the  influence  of  the  protu¬ 
berance  oersisting  downstream  for  a  hundred  protuberance  heights  or  more. 

1  * 

Sedney  has  described  the  effects  of  two-  and  three-dimensional  protu¬ 
berances  on  boundary-layer  flows.  In  the  two-dimensional  case  at  low 
speeds,  transition  generally  occurs  downstream  of  the  reattachment  point 
if  the  protuberance  height  k  is  considerably  less  than  the  undisturbed 
boundary  layer  thickness  6. 

This  reoort  describes  the  flow  field  near  a  square  two-dimensional 
protuberance  immersed  in  a  plane  Couette  flow.  Numerical  solutions  are 
based  on  a  finite-difference  analog  of  the  incomoressible  Navier-Stokes 
equations.  The  solution  technique  is  an  extension  of  a  method  developed 

p 

by  Gosman  et_  al_.  for  the  solution  of  the  steady  state  stream  function 
and  vorticity  formulation  of  the  Navier-Stokes  equations. 

The  object  of  this  work  is  to  describe  the  steady  flow  disturbances 
created  when  a  two-dimensional  protuberance  is  immersed  in  a  laminar  flat 
plate  boundary  layer.  It  should  be  noted  that  the  boundary  layer  equa¬ 
tions  are  not  valid  for  describing  the  flow  near  the  protuberance,  but 
would  most  likely  be  aoolicable  far  downstream  of  the  reattachment  ooint. 
In  this  study  a  nlane  Couette  flow,  an  exact  solution  of  the  Navier- 
Stokes  equations,  is  used  to  model  the  undisturbed  flow  on  the  flat 
01  ate.  We  exDect  the  qualitative  features  of  the  separation  phenomena 
induced  by  the  protuberance  in  the  Couette  flow  to  be  similar  to  those 
found  in  the  flat  Dlate  boundary  layer  case. 

Figure  1  shows  a  schematic  of  the  flow  field  near  a  square  orotu- 
berance  immersed  in  a  plane  Couette  flow.  The  protuberance  is  mounted 
on  the  stationary  bottom  plate.  The  upper  olate  moves  at  a  constant 
speed  UQ.  Far  upstream  of  the  protuberance,  shown  schematically  as 

position  1,  the  velocity  distribution  is  prescribed  by  a  linear  Couette 
flow  profile.  As  the  flow  nears  the  protuberance,  we  expect  it  to 
separate,  forming  a  small  upstream  recirculation  region. 

Downstream  of  the  protuberance,  at  position  2,  we  expect  separated 
flow  profiles  with  the  height  of  the  reversed  flow  region  no  larger  than 
the  protuberance  height  k.  Velocity  profiles  at  this  position  will 
exhibit  a  velocity  overshoot  near  the  moving  plate  (u>U_)  due  to  the 

flow  constriction  caused  by  the  protuberance.  Far  downstream  of  the 


* References  are  listed  on  page  38, 
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RECOVERING"  WAKE  PROFILES 

Figure  1.  Square  Protuberance  Immersed  in  Plane  Couette  Flow 


protuberance  and  the  reattachment  point,  at  position  3,  the  flow  still 
exhibits  a  perturbation  from  the  original  Couette  flow.  These  perturba¬ 
tions  decay  and  the  flow  returns  to  the  Couette-tyne  distribution 
farther  downstream. 

Numerical  results  are  presented  which  describe  the  flow  near  such 
protuberances  for  Re  =  1,  10,  25,  50,  100  and  200;  with  the  Reynolds 
number  Re  based  on  the  plate  velocity  UQ  and  the  orotuberan.ee  height  k. 

The  local  roughness  Reynolds  number,  Rek  =  Ukk/v,  is  one-half  of  each  of 

the  above  values;  where  Uk  is  the  mean  velocity  that  would  exist  in  the 

undisturbed  flow  at  the  protuberance  height  if  the  protuberance  was  not 
present.  The  numerical  results  are  illustrated  with  CALCOMP  plots  of 
streamlines  and  vorticity  contours  in  the  flow  field  around  the 
protuberance. 


II.  FORMULATION  OF  THE  PROBLEM 


The  mathematical  description  of  the  disturbed  flow  near  the  protu¬ 
berance  should  be  based  on  the  Navier-Stokes  equations  because  viscous 
separation  phenomena  are  important.  The  recirculation  region  downstream 
of  the  protuberance  and  the  position  of  flow  reattachment  along  the 
stationary  plate  influence  the  length  of  the  recovery  zone  downstream 
of  the  protuberance.  In  the  two-dimensional  case,  the  Navier-Stokes 
equations  can  be  expressed  in  rectangular  coordinates  by  the  x-momentum 
equation 


,,9u  ,  „3u  _  1  aP  t 

IF —  +  V -  =  - - +  \) 

ax  ay  P  ax 

and  the  y-momentum  equation 


3V  +  v9V._l!P+v( 

f  'Nv/  ~  'N  \/  v  v 
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The  continuity  equation  is 
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7? 
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ax  ay 
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The  pressure  can  be  eliminated  from  the  momentum  equations  by  differen¬ 
tiating  equation  (1)  with  respect  to  y,  equation  (2)  with  respect  to  x, 
and  subtracting  the  resulting  equations.  The  continuity  equation  (3) 
can  be  satisfied  by  defining  the  stream  function  ip  such  that 


u 


f  and 


V 


a^ 

ax  * 


By  introducing  the  definition  of  vorticity  c,  where 


(4) 
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9v_  au_  . 
c  "  ax  "  ay  * 


(5) 


equations  (1),  (2)  and  (3)  can  be  combined  with  (4)  and  (5)  into  two 
equations:  the  ell i otic  vorticity  transport  equation 


and  an  elliptic  Poisson-type  equation 

4+4+^  =  °-  (7) 

3x  ay 


Equations  (6)  and  (7)  represent  the  stream  function  and  vorticity 
formulation  of  the  Navier-Stokes  equations.  These  equations  govern  the 
flow  subject  to  the  prescribed  boundary  conditions.  It  should  be  noted 
that  only  the  prescribed  boundary  conditions  and  flow  parameters  make 
one  numerical  solution  different  from  another.  The  accurate  specifica¬ 
tion  of  boundary  conditions  is  thus  very  important.  Since  these 
governing  equations  are  elliptic,  boundary  conditions  must  be  soecified 
on  all  boundaries  of  the  computational  region. 


III.  BOUNDARY  CONDITIONS 


7 

Roache''  discusses  many  methods  for  treating  boundary  conditions  for 
the  stream  function  and  vorticity  in  numerical  solutions  of  the  Navier- 
Stokes  equations.  It  should  be  noted  that  many  of  the  methods  which 
have  been  used  in  the  past  to  prescribe  boundary  conditions  for  the 
vorticity  were  inconsistent  and/or  unsuccessful.  The  downstream  or  out¬ 
flow  boundary  and  the  share  protuberance  corners  have  been  particularly 
troublesome  in  many  instances.  The  difficulty  at  the  sharp  corners  will 
be  discussed  in  detail. 

Figure  2  shows  an  outline  of  the  computational  area  employed  in 
the  numerical  calculations  and  the  notation  used  along  boundaries.  A 
plane  Couette  flow  is  specified  along  the  upstream  boundary  AH  so  that 

u  =  U0(y/2k)  and  v  =  o;  (8) 

with 

*  and  c  =  Uo/2k  (9) 

The  no  slip  condition  (u  =  v  =  o)  is  specified  along  the  stationary 
plate,  AB  and  EF,  and  the  protuberance,  BC,  CD  and  DE.  The  stream 
function  y  is  specified  to  be  a  constant  on  these  boundaries;  taken 
equal  to  zero  for  convenience.  The  vorticity  c  on  these  boundaries  can¬ 
not  be  specified  independently  of  the  flow  off  the  surface.  Appendix 
A  shows  a  derivation  of  the  relationship  between  the  vorticity  at  the 
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¥  from  Eqn.O  I) 

5  from  Eqn.(A7)  or(A8) 
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Figure  2.  Schematic  of  Computational  Region  and  Boundary  Conditions 


wall  and  the  flow  adjacent  to  the  wall*  Two  different  expressions  for 
the  wall  vorticity  have  been  used  in  obtaining  numerical  solutions  to 
examine  the  influence  of  the  wall  vorticity  equation  on  the  calculated 
results.  These  expressions  are  given  by  equations  (A7)  and  (A8). 

Appendix  B  describes  two  alternatives  for  evaluating  the  vorticity 
at  the  sharp  corners,  C  and  D,  on  the  protuberance.  The  sharp  concave 
corners  B  and  E  are  no  problem  because  c  =  0  at  these  points.  In  the 
present  difference  procedure,  the  value  of  c  at  B  and  E  doesn't  enter 
the  iteration  formulae  for  nearby  points,  but  its  value  is  needed  for 
vorticity  contour  plots. 


Along  the  downstream  or  outflow  boundary,  FG,  the  conditions 


dip 

ax 


0  and 


(10) 


are  specified.  The  Couette  flow  solution  is  one  member  of  the  family 
of  profiles  that  satisfies  these  boundary  conditions.  A  comoutational 
experiment  is  needed  to  determine  the  extent  of  upstream  influence  of 
the  errors  caused  by  specifying  equations  (10)  at  a  fixed  position 
downstream  of  the  Drotuberance.  Results  from  such  tests  will  be 
discussed. 


The  upper  plate,  GH,  moves  at  a  constant  velocity  UQ*  The  no  slip 
conditions  (u  =  UQ,  v  =  0)  are  specified  along  boundary  GH.  The  stream 
function  ip  is  thus  specified  by 

*  =  UQk  (11) 

and  the  vorticity  c  is  determined  by  using  either  equation  (A7)  or  (A8). 

IV.  NUMERICAL  PROCEDURE 

2 

Gosman  et  aK  describe  the  comnlete  set  of  finite-difference  equa¬ 
tions  which  were  used  in  this  numerical  study.  The  finite-difference 
representations  of  equations  (6)  and  (7)  are  based  on  the  "tank  and  tube" 
method.  The  flow  field  is  divided  into  an  array  of  cells  called  tanks 
which  are  connected  to  adjacent  cells  by  tubes.  Mass  is  convected  from 
one  cell  to  another  along  the  tubes.  The  difference  equations  are  not 
restricted  to  evenly  spaced  nodes  and  an  uneven  mesh  spacing  between  grid 
points  can  be  easily  specified  by  prescribing  the  x-  and  y-coordinate  of 
each  grid  point.  An  upwind  difference  scheme  is  used  for  the  convective 
terms  in  the  governing  equations  with  information  advected  into  a  cell 
only  from  cells  upwind  of  it.  The  finite-difference  formulation  should 

4 

therefore  be  classified  as  first-order  accurate.  Gosman  has  analyzed 
the  difference  equations  and  has  concluded  that  convergence  is  likely 
in  all  circumstances  with  constant  prooerty  flow.  Divergence  of  the 
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iteration  procedure  has  occurred,  however,  under  certain  circumstances. 


4 


The  iteration  orocess  used  in  this  study  is  the  Gauss-Seidel  proce¬ 
dure.  The  complete  grid  is  scanned  node  by  node,  and  the  values  of  ^ 
and  ^  at  a  particular  node  are  calculated  using  neighboring  values  from 
the  previous  iteration.  Calculations  in  the  first  iteration  are  based 
on  initial  guesses  for  \p  and  c  at  all  grid  points.  The  scanning  process 
is  continued  until  the  grid  is  completely  scanned.  Special  formulae  are 
necessary  on  the  boundaries.  This  orocedure  is  then  repeated,  with  each 
complete  scan  equal  to  one  cycle  of  iteration,  until  the  criterion  for 
convergence  is  satisfied  or  the  calculation  is  otherwise  terminated.  If 
convergence  is  achieved,  the  velocity  distribution,  skin  friction  coeffi¬ 
cient,  etc.  are  calculated  and  a  plot  routine  is  executed. 


Convergence 


is  tested  by  the  calculation  of  the  fractional  changes 


<  x, 


Maximum  in  field 


and 


<  X; 


Maximum  in  field 


(12) 


03) 


where  \pQ  and  cQ  are  selected  reference  values  of  \p  and  c,  typically  0.5 

and  1.0,  respectively.  The  value  of  x  was  set  at  0.0001  for  all  calcu¬ 
lations.  Equations  (12)  and  (13)  are  applied  by  calculating  the 
fractional  change  of  \p  and  c  at  every  grid  point  in  the  field  for  each 
iteration.  When  the  absolute  value  of  the  fractional  change  is  less 
than  x  at  every  point,  convergence  is  examined  in  more  detail.  Care 
must  be  taken  to  examine  the  convergence  of  each  case  because,  with  the 
fractional  change  criteria,  a  slowly  converging  process  may  be  mistaken 
for  one  that  is  fully  converged.  Calculated  values  of  \p  and  c  at 
selected  points  in  the  field  are  examined  over  a  period  of  approximately 
ten  iterations  to  insure  that  convergence  has  been  achieved  if  equations 
(12)  and  (13)  are  satisfied.  Convergence  was  not  achieved  in  the 
calculations  for  cases  with  Re  >  300. 


V.  DESCRIPTION  OF  NUMERICAL  CALCULATIONS 


Numerical  solutions  have  been  obtained  for  Re  =  1 ,  10,  25,  50,  100 
and  200  which  describe  the  flow  around  a  square  protuberance  immersed  in, 
a  plane  Couette  flow.  The  iteration  procedure  used  in  this  study  did 
not  converge  for  cases  which  were  attempted  at  Re  =  300,  500  and  1000. 
Computational  results  will  be  illustrated  with  selected  velocity  vector 
plots,  vorticity  contours  and  streamlines  in  the  flow  field  drawn  on  an 
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off-line  CALCOMP  plotter.  A  contour  plotting  subroutine  developed  by 

Hartwig^  was  used  for  interpolation  and  contour  plotting  control.  Calcu¬ 
lations  were  performed  in  double  precision  on  the  BRLESC  I  and  II 
computers  at  the  Ballistic  Research  Laboratories. 


Two  different  grid  mesh  configurations  were  used  in  the  calculations. 
The  grid  for  Re  <  50  consisted  of  a  6  x  11  cell  protuberance  immersed  in 
a  51  x  21  cell  fTow  field.  The  mesh  aspect  ratio  e  =  ax/a y  varied  from 
2.0  near  the  protuberance  to  12.0  at  the  downstream  boundary.  Calcula¬ 
tions  for  Re  >_  100  were  oerformed  using  a  6  x  11  cell  protuberance 
immersed  in  a  81  x  21  cell  flow  field  with  3  varying  from  2.0  near  the 
protuberance  to  42.6  at  the  downstream  boundary.  The  inflow  boundary, 

AH  in  Figure  2,  was  located  at  approximately  X  =  -8.0,  where 

x  =  (x-xk)/k  (14) 


refers  to  distances  measured  from  the  position  of  the  protuberance 
leading  edge,  x^,  in  terms  of  the  protuberance  height  k. 


The  outflow  boundary,  FG,  was  located  at  apnroximately  X  =  9.5  in 
the  51  x  21  grid  and  at  approximately  X  =  93.5  in  the  81  x 21  grid.  The 
mesh  spacing  in  the  y-direction  was  fixed  at  Ay  =  k/10.  The  mesh 
spacing  in  the  x-direction  varied  from  ax  =  k/5  near  the  orotuberance  to 
Ax  =  ( 1 . 2 ) k  at  the  outflow  boundary  in  the  51  x  21  grid  and  to  ax  =  (4.3)k 
at  the  outflow  boundary  in  the  81  x  21  grid.  The  mesh  spacing  downstream 
of  the  protuberance  was  stretched  according  to 


X 


i+1 


xi 


+  k/5  + 


1 tan(iA/1Q) [  . 
2  tan(A)  * 


(15) 


where  A  =  1.0  and  i  =  1,2,  “10  in  the  51  x  21  grid,  and  A  =  0.05  and 
i  =  1,2, “40  in  the  81  x  21  grid. 


Figure  3  shows  a  velocity  vector  plot  of  the  flow  around  the  protu¬ 
berance  at  Re  =  1.  The  location  of  grid  Doints  in  the  51  x  21  grid  can 
be  visualized  in  this  figure--a  grid  point  is  located  at  the  tail  of  each 
velocity  vector.  The  magnitude  of  the  plate  velocity  UQ  is  shown  by  the 

length  of  the  vectors  above  the  moving  plate.  All  vectors  are  scaled 
with  respect  to  UQ.  This  plot  shows  that  at  very  low  Re  the  flow 

disturbances  caused  by  the  protuberance  are  confined  to  a  zone  within  a 
few  protuberance  heights  on  either  side  of  the  protuberance.  The  flow 
returns  to  an  exact  Couette  flow  distribution  by  X  =  9.0;  thus  flow 
disturbances  are  confined  to  a  zone  very  near  the  protuberance  at  Re  =  1 . 


Figure  4  shows  vorticity  contours  and  streamlines  in  the  protu¬ 
berance  flow  field  for  Re  =  1 .  The  flow  conditions  are  the  same  in 
Figures  3  and  4.  The  arrows  shown  above  the  moving  plate  in  the  stream¬ 
line  and  vorticity  contour  plots  are  not  plotted  to  the  same  scale  as  in 
the  velocity  vector  plots.  The  numerical  treatment  of  the  boundary 
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CORNER  VORTICITY  (DOUBLE- VALUED) 


conditions  for  each  case  is  indicated  below  the  contour  plots.  Figure  4 
represents  results  obtained  using  a  second-order  formulation  for  the  wall 
vorticity  given  by  equation  (A8),  treating  the  vorticity  at  the  sharp 
convex  corners  using  the  double-valued  method  described  in  Appendix  B. 

Figure  4  shows  that  at  Re  =  1 ,  the  flow  field  is  nearly  symmetric 
about  the  mid-plane  of  the  protuberance,  X  =  0.5.  We  expect  flow  field 
symmetry  in  the  Stokes  flow  limit  as  Re->0.  Eleven  streamlines  are 
plotted  in  the  non-recirculating  flow  region  in  all  cases  which  will  be 
discussed.  It  should  be  remembered  that  the  upper  plate  is  one  stream¬ 
line  and  the  lower  plate  and  protuberance  is  another.  The  difference  in 
stream  function  between  two  adjacent  streamlines  is  specified  so  that 

A4,  ,  ’‘'upper  plate^'  *lower  plate  .  (16) 

where  A\p  is  proportional  to  the  mass  flow  rate  per  unit  area  between  two 
adjacent  streamlines.  Figure  4  shows  that  the  mass  flow  rate  between 
the  bottom  plate  and  the  adjacent  streamline  is  much  lower  than  between 
any  other  two  adjacent  streamlines,  as  a  consequence  of  the  low  flow 
velocities  near  the  stationary  plate. 


Streamlines  in  the  separation  bubbles  upstream  and  downstream  of 
the  protuberance  have  negative  values  of  ip,  with  |i|;|  very  small.  A 
different  spacing  is  used  in  these  regions,  depending  on  the  strength 
of  the  recirculation.  Either  one  or  as  many  as  four  streamlines  are 
plotted  in  the  separation  bubbles  depending  on  the  magnitude  of  i p. 

For  the  case  shown  in  Figure  4,  the  downstream  separation  bubble 
consisted  of  three  grid  points  with  negative  values  of  stream  function. 
The  plotting  routine  denoted  this  with  a  single  closed  contour  as 
indicated  in  the  figure.  The  upstream  separation  bubble,  on  the  other 
hand,  was  much  smaller— consisting  of  a  single  grid  point  with  negative 
ip.  Since  the  apparent  size  of  the  separation  bubble  depends  on  the 
^-contour  plotted,  the  location  of  separation  and  reattachment  should 
be  determined  from  the  vorticity  contours  for  c  =  0. 


Vorticity  contours  in  the  flow  field  for  Re  -  1  are  also  shown  in 
Figure  4.  Negative  values  of  vorticity  are  plotted  with 

=  cmax.  negative  value  in  field  (17) 

between  adjacent  vorticity  contours.  The  closely  spaced  contours  near 
the  protuberance  convex  corners  indicate  large  vorticity  gradients  at 
these  points.  Positive  values  of  vorticity,  indicating  counter¬ 
clockwise  rotation  of  fluid  elements,  are  found  adjacent  to  the  moving 
plate.  Four  vorticity  contours  are  plotted  in  this  region  with 


AC 


Siiax.  positive  value  in  field 


(18) 
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between  adjacent  contours. 

As  can  be  shown,  points  of  zero  vorticity  along  the  stationary  plate 
and  protuberance  surface  correspond  to  points  of  flow  separation  and 
reattachment.  The  sharp  concave  corners  at  the  protuberance-plate 
intersection  are  also  points  with  vorticity  equal  zero.  The  situation 
on  the  moving  plate  is  different.  Points  along  the  moving  plate  with 
vorticity  zero  do  not  indicate  the  position  of  separation  or  reattachment 
as  can  be  shown  from  equation  (A7).  A  contour  for  c  =  0  is  plotted  near 
the  moving  plate,  separating  the  contours  for  c  positive  from  those  for 
c  negative.  Contours  connecting  points  of  zero  vorticity  are  plotted  in 
Figure  4  and  can  be  used  to  locate  points  of  separation  and  reattachment 
by  comparing  these  contours  with  the  separated  flow  streamlines.  The 
predicted  location  of  separation  and  reattachment  depends  on  the  method 
used  to  treat  the  vorticity  boundary  conditions. 

Calculations  for  Re  =  10  illustrate  the  influence  of  the  vorticity 
boundary  conditions  on  the  numerical  solution.  Figure  5  shows  stream¬ 
lines  and  vorticity  contours  for  Re  =  10  using  the  second-order  formula¬ 
tion  for  wall  vorticity  (A8)  and  the  double-valued  treatment  at  the 
convex  corners.  The  upstream  separation  bubble  is  seen  to  be  slightly 
smaller  than  at  Re  =  1 ,  with  the  downstream  separation  bubble  much 
larger.  Reattachment  of  the  downstream  separation  bubble  occurs  at 
approximately  X.  =  3.0. 

The  numerical  solution  depicted  in  Figure  5  shows  an  inconsistency 
near  the  separation  point  on  the  downstream  protuberance  face.  If 
separation  is  defined  to  occur  where  c  =  0,  then  separation  occurs  less 
than  one  cell  below  the  downstream  corner.  If  separation  is  defined  by 
extrapolating  the  dividing  streamline  back  to  the  downstream  protuberance 
face,  then  separation  occurs  approximately  one  and  one-half  cells  below 
the  corner.  This  inconsistency  is  not  present  in  the  cases  for  higher 
Reynolds  numbers  and  at  Re  =  10  it  occurs  only  when  the  double-valued 
vorticity  method  is  used  at  the  sharp  convex  corners.  It  thus  apDears 
that  for  the  case  shown  in  Figure  5  we  can  infer  only  an  approximate 
position  of  separation,  namely  about  (0.1 )k  below  the  downstream  corner. 

Figure  6  shows  that  the  results  for  Re  -  10  are  practically 
unchanged  if  the  first-order  formulation  for  wall  vorticity  (A7)  is 
employed  instead  of  the  second-order  formulation  used  in  Figure  5.  The 
inconsistency  in  the  separation  point  location  is  still  present  with 
separation  approximately  (0.1 )k  below  the  downstream  corner. 

The  positions  of  separation  and  reattachment  are  shifted  slightly 
for  Re  =  10  if  the  single-valued  method  for  corner  vorticity  is  used 
instead  of  the  double-valued  method.  The  single-valued  method  biases 
the  results  toward  separation  at  the  sharp  corner.  Figure  7  shows  that 
separation  occurs  within  the  first  cell  below  the  corner,  with  no 
inconsistency  using  and  c  to  locate  separation.  Reattachment  now 
occurs  farther  downstream  at  X  =  3.2. 
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REYNOLDS  NUMBER  =  10 
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CORNER  VORTICITY  (DOUBLE-VALUED) 
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CORNER  VORTICITY  (OOUBLE- VALUED) 
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WALL  VORTICITY  (SECOND-ORDER) 
CORNER  VORTICITY  (SINGLE-VALUED) 

Second-Order  Single-Valued  Corner  Vorticity 


Figure  8  shows  that  the  position  of  downstream  separation  changes 
slightly  if  the  first-order  formulation  for  wall  vorticity  (A7)  is  used 
instead  of  second-order  formulation  shown  in  Figure  7.  The  position  of 
reattachment  remains  at  X  =  3,2,  Comparison  of  Figures  7  and  8  shows 
that  the  prediction  of  the  upstream  separation  bubble  is  affected  slightly 
by  the  order  of  accuracy  used  in  calculating  the  wall  vorticity. 

It  is  interesting  to  note  that  in  snite  of  the  small  changes  caused 
by  the  four  different  boundary  conditions,  the  general  features  of  the 
predicted  flow  fields  are  similar.  The  upstream  separation  bubble  is 
very  small--less  than  (0.5)k  in  height  and  length.  The  downstream 
separation  bubble  has  a  height  no  greater  than  the  protuberance  height 
and  a  length  equal  to  approximately  two  protuberance  heights  at  Re  =  10. 

As  has  been  shown,  the  locations  of  separation  and  reattachment  cannot 
be  resolved  precisely  because  of  the  uncertainty  about  the  boundary 
condi tions . 

Figure  9  shows  streamlines  and  vorticity  contours  in  the  flow  field 
for  Re  =  25  based  on  the  second-order  formulation  for  wall  vorticity  and 
the  double-valued  technique  at  the  corners.  The  upstream  separation 
bubble  is  very  small --approximately  (0 . 3 ) k  high  based  on  the  c  =  0 
contour.  The  downstream  separation  bubble  is  almost  ( 4 . 1 5 ) k  long  at  this 
Reynolds  number.  It  appears  that  as  the  Reynolds  number  increases,  the 
predicted  location  of  downstream  separation  is  not  as  sensitive  to  the 
vorticity  boundary  conditions  as  at  Re  =  10.  There  is  no  inconsistency 
in  locating  the  position  of  downstream  separation  at  Re  =  25.  Separation 
occurs  less  than  one  cell  below  the  downstream  corner.  Calculations 
were  performed  for  Re  =  25  using  other  wall  vorticity  boundary  conditions 
but  the  results  differ  only  slightly  from  those  shown  in  Figure  9. 

Figure  10  shows  velocity  vectors  in  the  flow  field  for  Re  =  50.  As 
discussed  previously,  a  velocity  overshoot  occurs  adjacent  to  the  moving 
plate  due  to  the  flow  constriction.  U-components  as  large  as  ( 1 . 44 ) UQ 

are  present  in  the  flow  field  at  Re  =  50.  Figure  10  shows  that  appre¬ 
ciable  flow  disturbances  exist  at  the  downstream  boundary  located  at 
X  =  9.5.  It  is  expected  that  the  x-wise  decay  of  these  disturbances 
could  be  predicted  if  the  grid  points  were  continued  beyond  the  present 
downstream  boundary.  The  recovery  zone  may  extend,  however,  for  several 
hundred  protuberance  heights. 

Figure  11  shows  streamlines  and  vorticity  contours  in  the  flow 
field  for  Re  =  50.  These  plots  are  based  on  a  numerical  solution  using 
the  second-order  formulation  for  the  wall  vorticity  and  the  double¬ 
valued  vorticity  technique  at  the  convex  corners.  A  very  small 
separation  bubble  is  predicted  upstream  of  the  protuberance.  The  flow 
separates  less  than  one  cell  below  the  downstream  corner  and  reattaches 
at  X  =  8.2. 
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REYNOLDS  NUMBER  =  10 
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Figure  8.  Contour  Plots  for  Re  =  10  with  First-Order  Single-Valued  Corner  Vorticity 
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’-‘igure  10.  Flow  Field  Velocity  Vectors  for  Re  -  50  with  Second-Order  Double-Vain 
Corner  Vorticity 
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CORNER  VORTICITY  (DOUBLE  -  VALUED) 


A  computational  test  was  performed  to  determine  the  influence  of  the 
downstream  boundary  position  on  the  results  shown  in  Figure  11.  The 
downstream  boundary  in  Figure  11  is  at  X  =  9.5.  This  boundary  was  moved 
farther  downstream  to  approximately  X  =  14.0  and  the  calculation  for 
Re  =  50  was  repeated.  There  was  no  change  in  the  results  to  the  scale 
of  the  contour  plots  shown  in  Figure  11. 

A  computational  test  was  also  performed  to  determine  the  influence 
of  the  inflow  boundary  position  on  the  results  shown  in  previous  figures. 
The  inflow  boundary  was  placed  at  X  =  -4.0  and  calculations  for  Re  =  10, 

25  and  50  were  repeated.  The  results  near  the  protuberance  only  changed 
slightly,  but  it  was  evident  that  perturbations  caused  by  the  presence 
of  the  protuberance  were  affecting  the  calculations  for  grid  points 
adjacent  to  the  inflow  boundary;  i.e.  grid  points  next  to  the  inflow 
boundary  no  longer  showed  an  exact  Couette  flow  behavior.  With  the 
inflow  boundary  at  X  =  -8.0  the  exact  Couette  flow  solution  is  predicted 
for  the  entire  region  between  the  inflow  boundary  and  X  =  -4.0  for 
Re  =  50. 

Calculations  for  Re>50  were  performed  using  the  larger  81  x  21  grid 
to  accommodate  the  large  separated  flow  region  expected  at  higher 
Reynolds  numbers.  The  downstream  boundary  was  located  at  X  =  93.5. 

Figure  12  shows  streamlines  and  vorticity  contours  in  part  of  the  flow 
field  for  Re  =  100.  Separation  occurs  less  than  one  cell  below  the 
downstream  corner  with  reattachment  at  X  =  14.2.  The  calculations  pre¬ 
dict  that  the  flow  disturbances  introduced  by  the  protuberance  have  not 
fully  decayed  at  X  =  93.5.  Perturbations  of  the  x-comoonent  of  local 
velocity  as  large  as  0.44%  of  the  Couette  flow  value  at  the  same  position 
above  the  stationary  plate  are  present  at  X  =  93.5.  The  largest 
perturbations  are  found  adjacent  to  the  stationary  olate--one  cell  height 
off  the  plate. 

Figure  13  shows  streamlines  and  vorticity  contours  in  part  of  the 
flow  field  for  Re  =  200.  Separation  occurs  less  than  one  cell  below 
the  downstream  corner  with  reattachment  at  X  =  26.0.  At  Re  =  200  the 
flow  disturbances  introduced  by  the  protuberance  decay  more  slowly  than 
at  Re  =  100.  At  Re  =  200  perturbations  of  the  x-component  of  local 
velocity  as  large  as  8.1%  of  the  Couette  flow  value  at  the  same  y-position 
are  present  at  X  =  93.5.  The  calculation  for  Re  =  200  required  382 
iterations,  using  crude  initial  guesses,  to  satisfy  the  fractional  change 
criteria  given  by  equations  (12)  and  (13).  This  required  approximately 
31  minutes  of  computer  time  on  BRLESC.  Cases  for  lower  Reynolds 
numbers  required  less  computer  time  and  fewer  iterations  for  convergence 
using  the  same  initial  guesses. 

Cases  for  Re  =  300,  500  and  1000  were  run,  but  convergence  was  not 
achieved  for  any  of  these  cases  in  500  iterations.  In  these  cases  the 
predicted  values  of  stream  function  and  vorticity  at  a  given  point  in 
the  flow  began  to  fluctuate,  with  no  apparent  pattern  toward  convergence, 
after  approximately  400  iterations.  Each  case  was  run  for  500  iterations 
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Figure  12.  Streamlines  and  Vorticity  Contours  Near  Protuberance  for  Re 
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and  then  terminated  for  non-convergence.  The  reasons  for  non-convergence 
could  not  be  determined.  It  is  interesting  to  recall  the  experiences  of 

Macagno  and  Hung6  and  Mueller  and  O'Leary^  in  related  problems. 

Macagno  and  Hungw  studied  the  flow  of  a  liquid  through  an  axisym- 
metric  condui t  expansion  by  numerical  calculations  and  used  exoerimental 
results  to  verify  their  calculations.  Their  computational  simulation 
was  based  on  two  different  approaches.  The  first  involved  an  iterative 
numerical  scheme  for  the  steady  flow  equations,  similar  to  the  approach 
used  in  the  present  study.  Their  second  approach  was  based  on  the 
equations  for  unsteady  flow,  retaining  the  local  acceleration  terms. 

Their  results  show  the  steady-state  iterative  procedure  becomes 
"unstable"  as  the  Reynolds  number  is  increased.  We  discuss  their  results 
in  terms  of  a  Reynolds  number,  Rec,  based  on  the  step  height  and  the 

centerline  velocity  in  the  smaller  upstream  pipe.  They  succeeded  in 
using  the  iterative  steady  flow  approach  for  Rec<100  and  the  results  were 

in  agreement  with  results  based  on  the  unsteady  equations.  Their  results 
for  Rer  =  100  and  200  are  based  on  the  unsteady  aporoach  because  "... 

for  the  same  mesh  size,  the  steady  approach  would  have  been  unstable." 

Calculations  of  laminar  flow  over  a  two-dimensional  backsteD  by 

Mueller  and  O'Leary^  were  performed  using  an  explicit  finite-difference 
form  of  the  unsteady  equations.  Their  Reynolds  number  Re^  is  based  on 

the  frees tream  velocity  and  step  height.  They  comnare  their  numerical 
results  with  experimental  data  for  Ref.<100  and  present  numerical  results 

for  Reynolds  numbers  as  large  as  Re^  =  200.  They  experience  no  computa¬ 
tional  instabilities  in  these  cases.  Their  experimental  work,  however, 
shows  that  for  Ref>100,  "...  three-dimensional  effects  appear  just 

before  reattachment  in  the  form  of  helical  motion.  The  recirculation 
center  then  appears  to  become  increasingly  unstable  and  eventually 
appears  to  break  up  into  several  small  center  eddies.  The  flow, 
possibly  because  of  three-dimensional  effects  and/or  with  the  onset  of 
vortex  shedding,  is  increasingly  unstable  as  the  Reynolds  number  is 
raised. " 

The  convergence  problems  which  were  experienced  in  the  present 
study  for  Re>200  are  probably  related  to  the  stability  problems  dis¬ 
cussed  by  Macagno  and  Hung6.  The  grid  spacing  ax  in  our  81  x  21  grid 
becomes  very  large  downstream  of  the  protuberance  with  Ax/k  as  large 
as  4.26  at  X  =  93.5.  It  is  also  possible,  however,  that  our  convergence 
problems  are  related  to  a  physical  flow  unsteadiness  of  the  type 

described  by  Mueller  and  O'Leary^  in  their  experiments  for  Re^>100. 

This  intriguing  question  cannot  be  resolved  in  the  present  numerical 
study. 
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VI.  DISCUSSION  AND  CONCLUSIONS 


The  flow  around  a  square  protuberance  immersed  on  a  olane  Couette 
flow  has  been  studied  by  computational  simulation  for  Reynolds  numbers 
up  to  200.  The  numerical  solutions  predict  a  very  small  separated  flow 
region  upstream  of  the  protuberance  with  length  and  height  oractically 
independent  of  Re.  The  downstream  separation  bubble,  on  the  other  hand, 
has  a  height  of  order  k  and  a  length  which  increases  almost  linearly 
with  Reynolds  number. 

Figure  14  shows  that  the  position  of  reattachment  of  the  downstream 
separation  bubble  varies  almost  linearly  with  Reynolds  number  according 
to  the  numerical  calculations.  Two  sets  of  experimental  data  for 
backs tep-type  flows  are  shown  for  comparison.  Measurements  by  Macagno 

and  Hung6  for  an  axisyntnetric  conduit  exoansion  are  plotted  against 
Re.,  based  on  the  centerline  velocity  and  sten  height.  Measurements  by 
c  7 

Mueller  and  O'Leary  for  flow  over  a  two-dimensional  backstop  are  clotted 
in  terms  of  Ref,  based  on  the  frees tream  velocity  and  sten  height. 

'  o 

Klebanoff  and  Tidstrom  (not  shown  in  Figure  14)  found  reattachment  at  X 
between  30  and  40  for  a  cylindrical  rod  in  a  flat  plate  boundary  layer. 
Their  experimental  results  are  for  Reynolds  numbers  between  550  and  916 
based  on  protuberance  height  and  freestream  velocity.  Unfortunately,  it 
is  not  possible  to  compare  our  numerical  results  in  a  quantitative  sense 
with  any  of  these  experimental  measurements  because  of  the  difference  in 
geometries.  It  does  appear,  however,  that  the  separated  flow  region 
downstream  of  a  square  protuberance  may  be  significantly  longer  than,  that 
found  in  the  backsteo  type  geometry  at  comparable  Reynolds  numbers. 

The  present  results  predict  that  the  flow  disturbances  introduced 
by  the  protuberance  decay  more  slowly  as  the  Reynolds  number  increases. 
Significant  velocity  perturbations  are  found  at  X  =  93.5  in  our  calcula¬ 
tions  for  Re  =  200.  These  perturbations  tend  to  decrease  the  skin 
friction  along  the  plate  downstream  of  the  protuberance  from  the  value 
that  would  exist  in  the  undisturbed  flow. 

The  difficulties  in  treating  the  boundary  conditions  for  vorticity 
at  the  sharp  protuberance  corners  have  been  discussed.  Results  for  Re  = 
10  have  shown  that  it  is  not  possible  to  resolve  the  exact  location  of 
separation  on  the  downstream  protuberance  face.  Numerical  results  for 
Re>25  predict  separation  to  occur  less  than  one  cell  below  the  down¬ 
stream  corner.  Calculations  for  a  two-dimensional  backstep  geometry  by 

Mueller  and  O'Leary7  also  predicted  separation  below  the  corner  so  they 
investigated  the  location  of  separation  in  their  experiments;  they 
remark,  "...  although  suggested  in  some  photographs,  separation  below 
the  corner  could  not  clearly  be  seen."  It  thus  seems  appropriate  to 
interpret  our  numerical  prediction  of  separation  less  than  one  cell 
below  the  corner  for  Re>25  as  separation  at  the  corner  itself  within  the 
accuracy  of  the  finite-difference  calculation. 
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RE  ATTACHMENT 


Figure  14.  Variation  of  Downstream  Reattachment  Position 
with  Reynolds  Number 
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APPENDIX  A 


The  proper  boundary  conditions  for  viscous  flow  along  a  moving  or 
fixed  wall  are  the  “no  slip"  conditions.  If  vorticity  is  used  as  a 
dependent  variable  in  a  numerical  solution,  then  relations  between  the 
vorticity  at  points  along  the  wall  and  the  flow  adjacent  to  the  wall 
must  be  used  in  the  calculation  which  are  consistent  with  the  “no  slip11 
conditions.  More  than  one  equation  for  the  wall  vorticity  can  be 
derived  depending  on  the  truncation  error  that  can  be  tolerated. 


y 


X 


•  W+2 
•W+l 


Fiaure  A-l  Grid  Notation  for  Vorticity  Boundary  Conditions 
at  Wall 


Consider  a  wall  moving  at  constant  speed  UQ  shown  in  Figure  A-l. 

The  relations  which  will  be  derived  are  also  valid  for  a  fixed  wall 
when  UQ  is  set  equal  to  zero.  We  expand  the  stream  function  ip  to  ^w+1 

by  a  Taylor  series  expansion  out  from  the  wall  value  such  that 


Vi  •  *»  *  <«>§ I.  *  V  <*') 

w  ay  w  ay  w 

The  terms  in  equation  (AT)  can  be  rewritten  by  using  the  "no  slip" 
conditions 

uw  =  Uo  and  vw  =  °>  (A2) 

the  definition  of  stream  function 

» ■  If «»«»--  H  •  <“> 
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(A4) 


and  the  definition  of  vorticity 

av  au 
5  ax  ”  ay  * 

The  vorticity  at  the  wall  can  be  expressed  as  the  second  derivative 

?w =  -  Hi  (A5) 

ay  w 

by  combining  equations  (A2) ,  (A3)  and  (A4).  Equation  (AT)  can  be 
rewritten  as 

*«+i =  *w  +  (iy)uo  -  ^  cw  -  ^-3[— -V"~]  +  0(Ay4)-  (A6) 

Equation  (A6)  yields  different  equations  for  the  wall  vorticity 
deDending  on  the  order  of  truncation  of  the  Taylor  series.  The  lowest- 
order  accurate  equation  results  from  retaining  terms  up  to  order 
o 

(Ay)4*  in  equation  (A6)  and  solving  for  c 

w 

?w  =  "(iy?  C*w+1  ’  '  (Ay)U°]  +  0(Ay)'  (A7) 

Equation  (A7)  represents  the  first-order  formulation  for  the  wall 
vorticity.  The  second-order  formulation 

‘  [,»*1  -  *.  '  -  T1  *  ,(«2)  (A8) 

is  obtained  by  retaining  terms  up  to  order  (Ay)  in  equation  (A6). 

Expressions  for  the  wall  vorticity  c  can  be  derived  in  a  similar 

manner  having  higher-order  accuracy.  Since  the  finite-difference 
equations  used  for  interior  points  are  only  first-order  accurate, 
equation  (AS)  should  be  of  sufficient  accuracy.  It  should  be  noted 
that  when  cw+-j^w  equations  (A7)  and  (AS)  give  two  different  approxima¬ 
tions  to  the  wall  vorticity.  If  (Ay)  is  used  to  represent  the  grid 
spacing  normal  to  the  wall,  then  equations  (A7)  and  (AS)  remain  the  same 
regardless  of  the  wall  orientation  or  the  value  of  ip  on  the  wall. 
Computational  results  will  be  used  to  illustrate  the  effects  of  using 
these  two  different  expressions  for  wall  vorticity. 
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APPENDIX  B 


The  description  of  the  wall  vorticity  at  the  concave  protuberance 
corners,  B  and  E  in  Figure  B-l,  presents  no  difficulty  because  c  =  0  ? 

at  these  points.  The  specification  of  vorticity  boundary  conditions  at 
the  sharp  convex  protuberance  corners,  C  and  D,  poses  a  difficult  comou- 
tational  dilemma.  A  certain  type  of  non-uniqueness  occurs  at  these  two 
corners  because  there  are  several  alternatives  available  for  evaluating 
the  wall  vorticity  at  these  points. 


Figure  B-l  Notation  and  Geometry  for  Vorticity  at  Sharp 
Corners 


Points  C  and  D  can  both  be  described  in  the  same  manner,  so  for  sake 
of  clarity  let  us  confine  our  attention  to  the  downstream  corner  D.  If 
point  D  is  considered  part  of  the  vertical  wall  DE,  the  first-order 
formulation  for  wall  vorticity  given  by  equation  (A7)  becomes 


(Bl) 


If  point  D  is  considered  to  lie  on  the  horizontal  wall  CD,  the  first- 
order  equation  for  wall  vorticity  becomes 


CD 


wi* 


^AA  ‘  *D^- 


(B2) 


Obviously,  only  in  special  cases  would  these  two  equations  be  equivalent. 
A  similar  set  of  equations  results  by  using  equation,  (AS)  representing 
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the  second-order  formulation  for  wall  vorticity; 


s 


D 


ana 


CBB 
2  ’ 


CD 


CAA 
2  • 


(B3) 


(B4) 


3 

Roache  suggests  that  both  equations  be  used.  If  the  first-order 
formulation  for  wall  vorticity  is  being  used,  (Bl)  and  (B2)  would  be 
employed  in  the  following  manner.  When  the  grid  scanning  procedure 
reaches  poyit  AA  just  above  point  D,  then  c D  given  by  equation  (B2)  is 

used  in  the  difference  equations  for  AA.  When  is  required  in  the 

difference  equations  for  point  BB,  equation  (Bl)  is  used.  This  technique 
for  calculating  the  corner  vorticity  treats  c  as  a  double-valued  function 
at  the  geometric  singularity  D.  This  is  one  of  two  techniques  that  have 
been  employed  in  the  nresent  calculations.  We  refer  to  this  as  the 
double-valued  technique.  A  similar  procedure  is  employed  with  equations 
(B3)  and  (B4)  if  the  second-order  formulation  is  being  used  to  calculate 
wall  vorticity. 


There  are  many  other  possibilities  for  treating  the  vorticity  at  the 

9 

sharp  convex  corners.  Roache  and  Mueller  have  investigated  seven 
different  methods  for  treating  the  vorticity  at  the  corner  in  a  backs tep- 
type  geometry.  Their  conclusion  is  that  at  higher  Reynolds  numbers  the 
calculated  results  are  not  very  sensitive  to  the  method  used. 

A  second  method  was  tested  in  our  calculations  in  an  attempt  to 
force  separation  to  occur  at  the  sharp  downstream  corner  D.  This  approach 
treats  ^  as  a  single-valued  function,  assuming  that  the  flow  along  wall 

CD  determines  the  value  of  the  vorticity  at  the  corner.  cn  is  calculated 

using  either  equation  (B2)  or  (B4),  consistent  with  the  procedure  being 
used  for  the  other  wall  points.  This  method  will  be  referred  to  as  the 
single- valued  technique.  Point  C  was  treated  as  though  it  were  on  the 
downstream  wall  CD  whenever  the  single-valued  technique  was  used. 

Our  experience  with  these  two  methods  showed  that  the  results  are 
sensitive  to  the  technique  used  at  low  Reynolds  numbers*  The  single¬ 
valued  technique  was  not  successful  at  forcing  separation  to  occur  at 
the  downstream  corner  although  the  point  of  separation  is  very  close  to 
the  corner--less  than  one  cell  below  the  corner. 
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